#!/usr/bin/perl -w
use strict;
use Bio::DB::Fasta;

my $dbh = Bio::DB::Fasta->new('genomes');
my %data;
my %strains;
for my $dir ( qw(JEL423.JEL423 JEL423.JAM81 JAM81.JEL423 JAM81.JAM81) ) {
    my ($target,$base) = split(/\./,$dir);
    $strains{$target}++;

    opendir(DIR, $dir) || die $!;
    for my $file ( readdir(DIR) ) {
	next unless ( $file =~ /(\S+)\.gff3$/ );

	open(my $fh => "$dir/$file" ) || die $!;
	while(<$fh>) {
	    next if /^\#/;
	    chomp;
	    my ($seqid,$src,$type,
		$start,$stop, $score,
		undef,undef,
		$group) = split(/\t/,$_);
	    my %group = map { my ($key,@vals) = split(/=/,$_);
			      $key,\@vals } split /;/,$group;
	    
	    $data{$base}->{$seqid}->{$start}->{$target}->{alleles} = join("/",@{$group{'alleles'}});
	    $data{$base}->{$seqid}->{$start}->{$target}->{score} = $score;
	}
    }
}

my @strains = sort keys %strains;

while( my ($base,$chroms) = each %data) {
    open(my $fh => ">$base.SNPs_giga.dat") || die $!;
    print $fh join("\t", qw(CHROM POS REF), @strains),"\n";
    for my $chrom ( sort keys %$chroms ) {
	for my $pos ( sort { $a <=> $b } keys %{$chroms->{$chrom}} ) { 
	    my $ref_allele = $dbh->seq($chrom,$pos,$pos);
	    print $fh join("\t", 
			   $chrom, $pos, $ref_allele,
			   (map { $chroms->{$chrom}->{$pos}->{$_}->{alleles} ||  $ref_allele}
			   @strains)),"\n";
	}
    }
}
